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The density of electronic one-particle states in monolayer graphene is studied by performing the 
Hybrid Monte-Carlo simulations of the tight-binding model for electrons on the n orbitals of car- 
bon atoms which make up the graphene lattice. Density of states is approximated as a derivative 
of the number of particles over the chemical potential at sufficiently small temperature. Simula- 
tions are performed in the partially quenched approximation, in which virtual particles and holes 
have zero chemical potential. It is found that the Van Hove singularity becomes much sharper 
than in the free tight-binding model. Simulation results also suggest that the Fermi velocity in- 
creases with interaction strength up to the transition to the phase with spontaneously broken chiral 
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Figure 1: A cut of the dispersion relation (left) and the density of states (right) for the non-interacting tight- 
binding model with different values of the staggered potential m. The inset on the left plot shows the filling of 
the hexagonal Brillouin zone of graphene with discrete lattice momenta on the 18x18 lattice with periodic 
boundary conditions and the positions of the Dirac point (K), Van Hove singularity (M) and the edge of the 
valence band (T). These points are also marked on the right plot (for m — 0). Solid lines with points on 
them illustrate the convolution (1.2) of the density of states on the 18x18 lattice with the derivative of the 
Fermi factor at T /k = 0.28. Lines correspond to the results of exact calculation and points - to the numerical 
results obtained with the discretized fermionic operator (see below). 

1. Introduction 

Electron transport properties of graphene are of great importance both for the industrial appli- 
cations of this novel material as well as for our understanding of the physics of strongly correlated 
electrons. It turns out that at low energies electronic excitations in graphene effectively behave as 
free massless Dirac fermions [1]. However, they propagate with the Fermi velocity v^- = c/300, 
which is much less than the speed of light c. There is also a natural way to make these Dirac 
fermions massive: one should break the symmetry between the two simple rhombic sublattices of 
hexagonal graphene lattice by introducing the so-called "staggered potential", so that the potential 
is lifted by some value m on the sites of one sublattice and lowered by m on the sites of other 
sublattice. This "staggered potential" m then plays the role of the Dirac mass. At higher energies 
of order of the hopping energy fc « 2.8 eV the quasiparticle dispersion relation becomes nonlinear. 
At E = ±V K 2 + m 2 , it has a saddle point, which results in the logarithmic divergence in the density 
of states. This Van Hove singularity leads to the so-called Lifshitz transition as the Fermi energy 
of the system approaches the saddle point. The total width of the valence band is 2\/9k 2 + m 2 . 
The dispersion relation and the density of states at different values of the staggered potential m are 
illustrated on Fig. 1 . 

Charge carriers in graphene are also subject to electromagnetic interactions. The treatment of 
these interactions can be significantly simplified due to the smallness of the Fermi velocity. First, 
retardation effects can be neglected and, second, the Lorentz force between the electrons also is 
suppressed by a factor vp jc as compared to the electric force. Thus it is sufficient to consider only 
the instantaneous Coulomb interaction between Dirac quasiparticles. 

An important question is how the parameters of the free tight-binding model such as the widths 
of the valence band, the Fermi velocity and the mass gap m change when one takes into account 
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the interactions between quasiparticles. In order to pose such a question, one has to adapt the 
Fermi-Landau liquid model and assume that despite the interactions, in some sense quasiparticles 
can be still described as free fermions with some modified dispersion relation. Strictly speaking, 
graphene is an example of the so-called marginal Fermi liquid [2], for which the Fermi-Landau 
liquid approximation becomes inapplicable in the vicinity of the Fermi point (Fermi point is a 
degenerate Fermi surface with zero radius) . However, as discussed in [3], even a small doping 
brings monolayer graphene back into the Fermi-Landau liquid regime. One can expect that the 
introduction of the "staggered potential", which eliminates the Fermi point, should have the same 
effect. 

In this paper the density of one-particle states in monolayer graphene is studied numerically. 
It is found that the Lifshitz phase transition, associated with the Van Hove singularity around the 
saddle point in the dispersion relation of the tight-binding model, becomes much more pronounced 
in the interacting theory. Such behavior is in agreement with the analytical calculations based on 
renormalization-group arguments [4]. There are also indications that the Fermi velocity increases 
with the interaction strength up to the transition to the phase with spontaneously broken chiral 
symmetry, again in agreement with renormalization-group arguments [5]. On the other hand, it 
seems that the width of the artificially induced (with the help of the staggered potential) mass gap 
remains practically constant. However, measurements with a better energy resolution are required 
to clarify the situation completely. 

In analytical calculations, it is most convenient to study the quasiparticle dispersion relation 
and the corresponding density of states by analyzing the poles of the fermionic Green functions. In 
numerical Monte-Carlo simulations, however, fermionic Green functions can only be calculated for 
a finite number of discrete values of imaginary (Euclidean) time. Extraction of the real-time Green 
functions from such Euclidean field correlators is an ill-defined numerical problem, which has no 
unique solution. The most commonly used method for extracting real-time quantities from Eu- 
clidean correlators is the Maximal Entropy Method (MEM) [6]. However, this method introduces 
significant systematic uncertainties, in particular, due to ambiguity of the choice of the "model 
function" which incorporates our prior knowledge about the expected form of the corresponding 
spectral function. Typically, MEM tends also to smear the singularities (such as thresholds or Van 
Hove singularities) of the spectral functions at the scales of order of temperature. An attempt to 
extract the AC conductivity of graphene from numerically calculated current-current correlators 
was already reported by the author in [7]. It was found that MEM indeed was not able to reproduce 
the AC conductivity with precision which would be sufficient for quantitative analysis. 

It should be noticed, however, that in a strict sense the density of states which can be formally 
extracted from the fermionic Green function does not correspond to the density of any eigenstates 
of the interacting Hamiltonian. The reason is that the notion of quasiparticle is not strictly defined 
in this case. The interpretation of the spectral function of interacting fermionic gas in terms of 
the density of quasiparticle states only makes sense in the framework of the Landau-Fermi liquid 
model. However, in this framework one can also think of many other possible definitions of the 
density of states, which are only constrained by the requirement to correctly reproduce the density 
of states for a free Fermi gas with an arbitrary dispersion relation for fermions. 

In this work, the density of states is defined as the derivative of the number of particles over 
chemical potential. For a gas of free fermions with the density of one-particle states p (E), the 
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particle density n (jit) and its derivative can be written as: 

n(jl)= fdEp(E) ~JW\ (L1) 

dn ^l = ^[dEp(E) l —— (1.2) 

^ 4cosh 2 (^) 

In the limit of zero temperature, the Fermi factor F _„ . becomes the Heaviside step func- 

1+exp (^) 

tion 6 (jli — £), and its derivative with respect to the chemical potential becomes the 5-function 
8{jl—E). Thus in the limit of zero temperature the derivative dn (ju) jd\i is exactly the density 
of one-particle states. At finite temperature the 5-function is smeared over a finite energy range 
with the width of order of temperature and becomes an exponentially decaying function. A direct 
estimate of the distribution width yields: 



-!- f dE — ^ p-^- = " v " ' (1.3) 

4cosh 2 (^ X 



(E-n_Y_ _n 2 (kTY 

E-n 
2kT 

Thus by performing simulations at sufficiently small temperatures and by measuring the derivative 
dn J^ , one can obtain a reasonably accurate estimate of the density of one-particle states. The 
resulting function is smeared as compared to the zero-temperature limit, similarly to the results 
which can be obtained by using MEM. On the other hand, such an estimate is not based on any 
"model function" and is thus free from possible bias in the measurements caused by the particular 
choice of this function. 

In practice, a direct numerical implementation of such measurements is a formidably difficult 
task due to the so-called "sign problem" of the Monte-Carlo simulations at finite chemical potential 
[8]. At zero chemical potential (that is, for graphene at half-filling) the sign problem is absent due to 
the symmetry between particles and holes [7, 9, 10]. For this reason, in this work the quantity ^jp- 
is calculated in the partially quenched approximation, in which virtual particles are not influenced 
by finite chemical potential. Such approximation can be introduced as follows: consider the tight- 
binding model with 1 + Nf independent fermion flavors and assume that chemical potential is 
nonzero only for Nf fermions. For the time being it is convenient to assume that the staggered 
potential m is equal to zero. Transforming the partition function of such a system into the path 
integral representation along the lines of [7, 9, 10], it is straightforward to obtain the following 
result: 

3T(n,T,N f ) = J | det (Af) 1 2 det (M + p ) Nf det (M - ju f f e~ s[ ^ , (1.4) 

where <p = (j)(X,z) is the Hubbard-Stratonovich [10] or the electrostatic potential [7, 9] field 
with the action S[(j>], x £ 0, (kT)^ 1 is the Euclidean time, X variable labels the sites of the 
graphene hexagonal lattice, M = d x — ho + i<f> (X, x) is the fermionic hopping operator and ho is 
the one-particle Hamiltonian of the tight-binding model. Consider now the linear response of 
the system to adding more flavours of fermions at finite chemical potential, that is, the derivative 
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dl0g ^dNf rN> ^ Wf-tO- This yields the partially quenched partition function 

2f q (H,T) = 3f- 1 (H = 0,T,N f = 0) y"^0|det(M)| 2 2ReTrlog(M + i u)^ SW . (1.5) 

This expression was derived by taking into account the invariance of the path integral under the 
reflection <p — > — <p and the anti-unitary symmetry of the operator M at zero chemical potential: 
CMC = — M 1 " with C = C(X,T\\Y, T2) = ±<5xy<5 (fi,^), where the plus and the minus signs are 
taken for even and odd sites of the graphene hexagonal lattice, respectively. Clearly, the square of 
C is just the identity operator. At low energies, when quasiparticles in graphene can be described 
as Dirac fermions, C plays the role of the 75 Dirac matrix. 

Finally, taking the derivatives of the quenched partition function over the chemical potential, 
one finds the partially quenched approximation for n (ju) and dn{\i) j d\i: 

n q (n) = 3f- 1 (n = 0,T,N f = 0) J #0|det(M)| 2 ^ReTr(M + /i) _1 e~ s ^ (1.6) 

= i r- 1 ( M =0,T,N f = 0) |^|det(M)| 2 ^ReTr(M + i u)" 2 e- s W, (1.7) 

where V is the number of lattice sites in the system. These expressions differ from the expression 
for the full, unquenched functions n(fi) and dn(jx) /dfX by the absence of chemical potential in 
the factor |det(M) | 2 . However, since the density of one-particle states is not a rigorously defined 
quantity in the interacting theory, there is no reason to believe that dn (ju) /dpL\^ = E is a better esti- 
mate of p (E) then dn q (fi) /d/J,\n = g. Indeed, for a free fermion gas both definitions are completely 
equivalent. 

The observables (1.6) and (1.7) are now well-suited for numerical simulations and can be cal- 
culated by averaging the quantities (M + jj, ) ~~ 1 and (M + /I ) ~ 2 over the equilibrium ensemble of the 
fields (X,t) with the weight |det(M) \ 2 e ~ s ^ [7, 9, 10]. To this end Hybrid Monte-Carlo simula- 
tions of the tight-binding model of graphene on the lattice with toric topology which consisted of 
18x18 hexagonal cells were performed. Temporal size of the lattice is L T = 18 with lattice spacing 
kAt = 0.2. Coulomb interactions are modelled by coupling the tight-binding model to the U (1) 
noncompact gauge field in 3 + 1 dimensions, as in [7]. The corresponding coupling constant is 
controlled by the dielectric permittivity £ of substrate on which the graphene monolayer is placed: 
a = 2ao/ (vf (e + 1)), where Ofo = e 1 /An = 1/137 [7]. Lattice size in the direction perpendicular 
to the graphene plane is L, = 18. Simulation setup and the discretization of the fermionic operator 
M are the same as in [7]. Despite the fact that the estimate dn '^ of the density of states in (1.7) 
can be expressed in terms of the partially quenched partition function (1.5) only at zero staggered 
potential m, the simulations have to be performed at some nonzero m (in this work m/K = 0.1 
and m/jc = 0.5) in order to ensure the invertibility of the fermionic operator M. The derivative 
dn fj? is found by first calculating n q (ju) according to (1.6) for equidistant values of /I separated 
by A/1 = 0.2 K and then numerically differentiating this function using the symmetric difference 

dnq(li) _ n q (ji + Aii)-n q (ii-An) 
d}X ~ 2A^i 

The resulting dependence of dn q (n) jd\i on \i is illustrated on Fig. 2 for ra/fc = 0.1 and 
m/K = 0.5. For comparison, the function dn q jd\l for the free theory and the density of states 
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Figure 2: Estimate (1.7) of the density of states for the 18 x 18 lattice with kAt = 0.2 (T/k = 0.28) and 
different values of the substrate dielectric permittivity at m/K = 0.1 (above) and m/K = 0.5 (below). Solid 
lines connecting the data points are cubic splines which are plotted to guide the eye. 



PDimc (E) for free Dirac fermions (which corresponds to dn q (ju) jd\i in the limit of zero tempera- 
ture) are also plotted: 

PDimc [E) = — — y \E\, \E\>m, p Dirac (E) =0, \E \< m. (1.9) 

This expression takes into account that on the hexagonal lattice with periodic boundary conditions 
the number of states in the element d 2 k of the momentum space is given by dN/V = d 2 k [7]. 

First of all one can note the sharp rise of the peaks at E w ±k, which are associated with the 
Van Hove singularities in the free theory. This is an indication that the Lifshitz phase transition 
(associated with the crossing of the Van Hove singularity by the Fermi level) might become much 
sharper in the interacting theory. Indeed, such phase transitions are known to be unstable with re- 
spect to even small interaction between particles. Application of renormalization-group techniques 
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to quasiparticles in the vicinity of the Van-Hove singularities also show that as the interactions are 
switched on, the density of states tends to increase [4]. In other words, the quasiparticle dispersion 
relation becomes more flat in the vicinity of the M point. This sharpening of the Van Hove singu- 
larity can be clearly seen for both values of m starting from the smallest nonzero coupling constant 
(which corresponds to substrate dielectric permittivity e = 10.0). It is also interesting to note a 
small asymmetry between the Van Hove singularities at E = — K and at E = K and a slight shift of 
the position of the minimum of the density of states from E = to positive values of E. 
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Figure 3: Minimal value of the estimate d "^^ of the density of states for ji in the range — K < ji < K for 
different values of the staggered potential. 

Another important prediction of the renormalization-group analysis, which has been recently 
confirmed experimentally [11], is the logarithmic divergence of the Fermi velocity in the interacting 
tight-binding model [5]. As follows from (1.9), the increase of Fermi velocity should result in a 
depletion of the density of states near E = 0. Superficially, since the slope of dn q (fi) jd\l near 
\l = clearly increases towards smaller e, it seems that our results point to the decrease of the 
Fermi velocity with interaction strength. However it might well be that due to the smearing (1.2) 
this increase of the density of states is just caused by the strong growth of the density of states 
at the Van Hove singularity. In order to demonstrate that the density of states near E = indeed 
decreases, on Fig. 3 the minimal value of dn q (li) jd\i for fi in the range — K < pL < K is plotted 
as a function of e. This minimum is always situated close to E = 0. One can see that indeed 
this minimal value decreases as the Coulomb interaction becomes stronger for e > 4. This is an 
indication of the decrease of the density of states in the vicinity of E = 0. As demonstrated in 
[7, 9], at smaller e the symmetry between simple sublattices of graphene hexagonal lattice breaks 
down spontaneously. In this phase there is no reason to believe that the Fermi-Landau liquid model 
is a good approximation. A glance at Fig. 2 suggest also that the width of the energy gap which 
is artificially induced by the staggered potential m is practically independent of the interaction 
strength. 

To conclude, the presented results are consistent both with the sharpening of the Van Hove 
singularity and with the increase of the Fermi velocity predicted by renormalization-group calcula- 
tions [4, 5]. It seems, however, that the former modification the dispersion relation is much stronger 
than the latter, and that much better energy resolution is required in order to study quantitatively 
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the evolution of the Fermi velocity and the mass gap in the tight-binding model with interactions. 
At present larger-scale simulations with significantly smaller temperatures (so that the smearing 
factor in (1.2) is much narrower) and with larger lattices are in progress, which would hopefully 
help to study the quasiparticle dispersion relation in more details. 
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